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We study lattice QCD at non-vanishing chemical potential using the complex Langevin equation. 

We compare the results with multi-parameter reweighting both from p. — 0 and phase quenched 
ensembles. We find a good agreement for lattice spacings below «0.15 fm. On coarser lattices the 
complex Langevin approach breaks down. Four flavors of staggered fermions are used on W = 4, 6 
and 8 lattices. For one ensemble we also use two flavors to investigate the effects of rooting. 
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I. INTRODUCTION AND OVERVIEW 


Dense and/or high temperature phases of strongly interacting matter are becoming experimentally accessible nowa¬ 
days due to heavy ion collision experiments at the Relativistic Heavy Ion Collider, the Large Hadron Collider, and 
especially the FAIR facility at GSI, as well as astrophysical observations of neutron starts. Theoretical understanding 
of the dense, strongly interacting phases and the first principles determination of the phase diagram of QCD as a 
function of the temperature and chemical potential are still lacking. This is a consequence of the sign problem, which 
makes lattice calculations at nonzero baryon density challenging. 

The standard non-perturbative tool for QCD, lattice QCD is defined by the path integral 


Z = J DUe-^^^detM{fj.) 


( 1 ) 


with the Yang-Mills action Sym of the gluons and the fermion determinant M{p) on a cubic space-time lattice. At 
nonzero chemical potential the determinant is non-real, therefore importance sampling methods are not applicable. 
For a review of ideas to circumvent the sign problem see EMI- 

One of the ways to avoid the sign problem is usingthe analyticity of the action and complexifing the field manifolds 
of the theory with the complex Langevin equation [J, (See also the related but distinct approach of the Lefschetz 
thimbles, where the integration contours are pushed into the complex plane @.) 

After promising initial results, it was noticed that the complex Langevin equation can also deliver convergent but 
wrong results in some cases [A 0] . Also technical problems could arise which are avoided using adaptive step-sizes 
for the L ang evin equation [3,1^ • In the last decade the method has enjoyed increasing attention related to real-time 
systems E3 t1J]i finite-density problems The method showed remarkable success in the case of 


finite densiW Bose gas [16| or the SU(3) spin-model 3^13, but the breakdown of the method was also observed a 
few times ESiHS]- The theoretical understanding of the successes and the failures of the method has improved: it 
has been proved that provided a few requirements (some ’offline’ such as the holomorphicity of the action and the 
observables, and some ’online’ such as the quick decay of the field distributions at infinity) the method will provide 
correct results [H: • 

It has been recently demonstrated that complex Langevin simulations of gauge theories are made feasible using 
the procedure of gauge cooling [13) HH ) see also [11] , which helps to reduce the fluctuations corresponding to the 
complexified gauge freedom of the theory. This method was first used to solve HDQCD (heavy dense QCD) where the 
quarks are kept static (their spatial hopping terms are dropped) [Hill, and it has been also extended to full QCD 
using light quarks in the staggered [Hj as well as the Wilson formulation [1^ . Gauge cooling makes the investigation 
of QCD with a theta term also possible [Hj. 

In this paper we compare results of the reweighting approach and the Complex Langevin approach for Np = 4 and 
Np = 2 QCD using staggered fermions. 

In Section|TTl we give a brief overview of the complex Langevin method. In Section lllll we summarize the reweighting 
method. In section lTVj we present our numerical results comparing the reweighting and complex Langevin simulations. 
Finally, we conclude in Section [V] 
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II. THE COMPLEX LANGEVIN EQUATION 


The Complex Langevin equation ii is a straightforward generalization of the real Langevin equation [s^. For 
the link variables of lattice QCD an update with Langevin timestep e reads [s^ : 


Ux,u{'^ + e) = exp 


^ ^ ^ “t” 's/^Vaxi') 


Ux,u{t), 


( 2 ) 


with Xa the generators of the gauge group, i.e. the Gell-Mann matrices, and the Gaussian noise rjaxi^- The drift force 
Kaxv is determined from the action 5'[f/] by 


Kax. = -Dax.S[U] (3) 

with the left derivative 

Dax.f{U)= (4) 

In case the drift term is non-real the manifold of the link variables is complexified to SL(3,C). The original theory is 
recovered by taking averages of the observables analytically continued to the complexified manifold. 

For the case of QCD the action of the theory involves the fermionic determinant through the complex logarithm 
function 


Seff = Sym — In detM{fj,). 


(5) 


The drift term in turn is given by 

Nf 

Kaxu = -Dax.SYM[U] + ^Tr[M-\^i,U)Dax.M{^i,U)\, ( 6 ) 

where the second term is calculated using one CG inversion per update using noise vectors [s^. The action we 
are interested in is thus non holomorphic, and in turn this results in a drift term which has singularities where the 
fermionic measure detM(^, f/) is vanishing. 

The theoretical understanding of the behavior of the theory with a meromorphic drift term is still lacking, but we 
have some observations as detailed below. Such a drift term seems to lead to incorrect results in toy models if the 
trajectories encircle the origin frequently [H, [s^ . In other cases the simulations yield a correct result in spite of a 
logarithm in the action [l^ |2C|. In an explicit example is presented where the simulations give correct results in 
spite of the frequent rotations of the phase of the measure. The condition for correctness is that the distribution of 
configurations vanishes sufficiently fast (faster than linearly) near the pole. 

For QGD itself we have a few indications that the poles do not affect the simulations at high temperatures: observing 
the spectrum of the Dirac operator [s^, comparisons with expansions which use a holomorphic action [^ . and the 
results presented in this paper. It remains to see whether simulations in the confined phase are affected. 

The ’distance’ of a configuration from the original SU(3) manifold can be monitored with the unitarity norm 

(7) 

X,fl 

where D = N^Nt is the volume of the lattice. In naive complex Langevin simulations, this distance grows exponentially, 
and the simulation breaks down because of numerical problems if it gets too large. This behavior can be countered 
with gauge cooling, which means that several gauge transformations of the enlarged manifold are performed in the 
direction of the steepest descent of the unitarity norm © 0,1111. With gauge cooling, the unitarity norm remains 
bounded at a safe level as long as the /3 parameter of the action is not too small. The value /3mm corresponds to a 
maximal lattice spacing, which seems to depend weakly on the lattice size, as can be checked easily for the cheaper 
HDQCD theory 0. 


III. REWEIGHTING 


In the multi-parameter reweighting approach one rewrites the partition function as [4 
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Z = J detM{fi) = j det M(/^o) t^M{^) } ’ 


where ^o is chosen such that the second line contains a positive definite measure which can be used to generate 
the configurations and the terms in the curly bracket in the last line are taken into account as an observable. The 
expectation value of any observable can be then written in the form: 


< O >0,^= 


Y,OW,f^)w{l3,l3o,fj.,fio) 

X;'u;(/3,/3o,/i,Aio) 


(9) 


with w{l3, Po, fj,, Ho) being the weights of the configurations defined by the curly bracket of eqn. ([5]). Note that 
gauge observables do not explicitly depend on /x, therefore their fi dependence comes entirely from the weight factors. 
Fermionic observables, on the other hand also explicitly depend on the chemical potential. 

In this paper we use two choices for the original, positive measure ensemble. The first choice is to use ho = Oj i-C- 
reweighting from zero chemical potential. For any choice of the target /3, /i parameters one can find the optimal /3o 
for which the fluctuation of the weights w{f3, /i) is minimal. This corresponds to the best reweighting line as discussed 
in We generated configurations at /r = 0 for /3 in the range 4.9 — 5.5. These were then used to reach the 

entire /3 plane via multi-parameter reweighting. Our second choice is to use the phase quenched ensemble, i.e. 
replacing det M(^o) by | det M(/i)| in eqs. dH]) and (ITUl) . In this case the reweighting factor contains only the phase of 
the determinant. 

For staggered fermions an additional rooting is required, for Np flavors the weights become 


w{P, /3o, /i, Ho) = 


det M (h) 


_detM{Ho) 


Nf/4 


( 10 ) 


Since for Np < 4 a fractional power is taken which has cuts on the complex plane it is important to choose these cuts 
such that the weights are analytic for real h values. This can be achieved by expressing det M{h) analytically as a 
function of h as discussed in [^, |4^ . 


IV. RESULTS 


We use the Wilson plaquette action for the gauge sector of the theory and unimproved staggered fermions with 
Np = 4 flavors if not otherwise noted. We have used three different lattice sizes for this study, 8^ x 4, 12^ x 6 and 
16^ X 8, all having the aspect ratio Ls/Lt = 2. 

Our main observables are the plaquette averages, the spatial average of the trace of the Polyakov loop 

J2TtP{x)/n!, P{x)= n U4ix,i) (II) 

X i=l..Nt 


and its inverse TrP ^{x)/Ng , the chiral condensate (ipip) the fermionic density n defined as 

1 /ainZ' 


(V'^) = 


I /ainZ 

\ dn 


( 12 ) 


fl \ dm 

with 17 the volume of the space-time lattice. We are also interested in the average phase of the fermion determinant, 
which measures the severity of the sign problem 

det M (h) 


(g2*c^) = 


det M(—/r) 


(13) 


We perform the complex Langevin simulations using adaptive step-size, with a control parameter which puts the 
typical step-sizes in the range e ~ 10“® — 5 x 10“^. Using such small step sizes allows us to avoid having to take the 
e —> 0 limit as the results are in the zero Langevin step limit within errors. We use initial conditions on the SU (3) 
manifold, and allow r = 10 — 30 Langevin time for thermalization, after which we perform the measurements for an 
other r = 10 — 30 Langevin time. We checked that proper thermalization is reached by observing that halving the 
thermalization time leads to consistent results. 

We have determined the pion masses as well as the lattice spacing using the wq scale as proposed in for several 
quark masses, see in Table ID One sees that choosing the quark masses ma = 0.05 for the Nt = 4 lattice, ma = 0.02 
for the Nt = 6 lattice and ma = 0.01 for the Nt = 8 lattice, in the vicinity of the critical temperature we have 
m^^/Tc ~ 2.2 — 2.4. We have additionally investigated the Nt = 8 lattice with am = 0.05, which corresponds to the 
rather heavy pion mass of m^^/Tc ~ 4.8. 
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/3 

arriq 

am-j^ 

a (fm) 

4.80 

0.01 

0.2458 ± 0.0007 

0.3355 ± 0.0001 

4.85 

0.01 

0.2480 ± 0.0007 

0.3315 ±0.0001 

4.90 

0.01 

0.2506 ± 0.0009 

0.3258 ± 0.0003 

4.95 

0.01 

0.2533 ± 0.0009 

0.3174 ± 0.0003 

5.00 

0.01 

0.2596 ± 0.0005 

0.2892 ± 0.0001 

5.05 

0.01 

0.2679 ± 0.0008 

0.2773 ± 0.0006 

5.10 

0.01 

0.2870 ± 0.0010 

0.1890 ± 0.0006 

5.15 

0.01 

0.3014 ±0.0012 

0.1410 ±0.001 

5.20 

0.01 

0.2957 ±0.0018 

0.1123 ±0.0006 

5.25 

0.01 

0.2918 ±0.0021 

0.0957 ± 0.0006 

5.40 

0.01 

0.2456 ± 0.0012 

0.0652 ± 0.0006 

4.80 

0.02 

0.3455 ± 0.0006 

0.3362 ± 0.0001 

4.85 

0.02 

0.3491 ± 0.0005 

0.3323 ± 0.0001 

4.90 

0.02 

0.3520 ± 0.0007 

0.3270 ± 0.0001 

4.95 

0.02 

0.3568 ± 0.0006 

0.3197 ±0.0001 

5.00 

0.02 

0.3629 ± 0.0007 

0.3084 ± 0.0002 

5.05 

0.02 

0.3726 ± 0.0005 

0.2878 ± 0.0004 

5.10 

0.02 

0.3898 ± 0.0008 

0.2420 ± 0.0009 

5.15 

0.02 

0.4073 ± 0.0006 

0.1751 ±0.0008 

5.20 

0.02 

0.4125 ±0.0011 

0.1341 ±0.0012 

5.00 

0.03 

0.4392 ± 0.0002 

0.2986 ± 0.0001 

5.10 

0.03 

0.4635 ± 0.0004 

0.2419 ± 0.0002 

5.20 

0.03 

0.4905 ± 0.0008 

0.1467 ± 0.0004 

5.25 

0.03 

0.4862 ± 0.0007 

0.1237 ± 0.0003 

5.40 

0.03 

0.4407 ± 0.0009 

0.0849 ± 0.0004 

4.80 

0.05 

0.5413 ± 0.0003 

0.3379 ± 0.0001 

4.85 

0.05 

0.5442 ± 0.0005 

0.3345 ± 0.0001 

4.90 

0.05 

0.5480 ± 0.0004 

0.3301 ± 0.0001 

4.95 

0.05 

0.5530 ± 0.0005 

0.3241 ± 0.0001 

5.00 

0.05 

0.5588 ± 0.0002 

0.3045 ± 0.0001 

5.05 

0.05 

0.5678 ± 0.0004 

0.3042 ± 0.0002 

5.10 

0.05 

0.5784 ± 0.0003 

0.2664 ± 0.0002 

5.15 

0.05 

0.5961 ± 0.0006 

0.2426 ± 0.0008 

5.20 

0.05 

0.6100 ± 0.0005 

0.1783 ± 0.0003 

5.25 

0.05 

0.6144 ± 0.0006 

0.1465 ± 0.0004 


TABLE I: Pion masses and lattice spacings for different /3 values and bare quark masses, measured on 12® x 24, 16® x 32 and 
24® X 48 lattices with Nf — 4. Statistical errors are indicated. 


A. Reweighting from fj, = 0 

First we have tested the theory at a fixed (3 = 5.4 at IVt = 4 as a function of /r, which is well above the deconfinement 
transition which at /r = 0 and m = 0.05 is at /3c ~ 5.04. 

In Fig. [1] we show the comparison of the gauge observables: plaquette averages and Polyakov loops. We generated 
O(IO^) independent configurations in the /r = 0 ensemble with the usual HMC algorithm (using every 50th configura¬ 
tion of the Markov chain), and we calculated the reweighting as detailed in section Hill One notes that the reweighting 
performs well for small chemical potentials n/T < 1 — 1.5, where there is a nice agreement between reweighting and 
OLE. The errors of the reweighting approach start to grow large as one increases fi above 1.5T, where the average of 
the reweighting is dominated by a few configurations. This is the manifestation of the overlap problem: the ensemble 
we have sampled has typical configurations which are not the typical configurations of the ensemble we wish to study. 

Next we turn to the fermionic observables: chiral condensate, fermionic density in Fig. [21 One notes that the 
reweighting of these quantities is possible to much higher values of /i/T. This is the consequence of their explicit 
dependence on /x, which dominates their change as the chemical potential is changed. This is in contrast to the gauge 
observables in Fig. [1] where the change is given entirely by the change in the measure of the path integral. The 
downward turn of the Polyakov loop and its inverse around /x/T = 3 is the result of the phenomenon of saturation: 
at this chemical potential half of all of the available fermionic states on the lattice are filled, as visible on Fig. |21 This 
lattice artifact can also be observed with static quarks [s^, and even in the strong coupling expansion [4^. 

Finally, the average phase factor in Fig.|2]is a good indicator of the severeness of the sign problem in the theory. One 
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FIG. 1: Comparison of plaquette averages and Polyakov loops and inverse Polyakov loops (defined in and below eg. (1111) 1 
calculated with OLE and reweighting from the /r = 0 ensemble. 




FIG. 2: Comparison of the chiral condensate and fermionic density as well as the phase average m calculated with CLE and 
reweighting from the = 0 ensemble. 



relative weight 

FIG. 3: The histogram of the relative weight of the configurations for different chemical potentials. 
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FIG. 4: Comparison of plaquette averages and Polyakov loops calculated with CLE and reweighting from the jj. — 0 ensemble. 
We use Nf = 2 fermion flavours for the data presented in this plot. 


sees that the average phase in the region /i/T> 1.5 — 2 indeed gets very small. Note that to see agreement between 
CLE and reweighting one has to be careful to choose the observable to be the analytic continuation of an observable 
on the SU(3) manifold. In this case we define the phase factor from the analytic continuation of the determinants, as 
written in (EJ. 

In Fig. [3] we show the histogram of the absolute value of the weights of the conhgurations normalized by the 
biggest weight in the ensemble. This illustrates the overlap problem: the ’further’ one tries to reweight from the 
original ensemble, the less and less will be the contribution of an average configuration to the average, which becomes 
dominated by very few configurations. Thus the fluctuations of the result become larger, and even the errorbars are 
not reliable as the distribution of the observables becomes non-Gaussian. As we show below, this situation improves if 
one chooses an ensemble ’closer’ to the target ensemble: in this case taking the phasequenched ensemble (|detM(/i)|) 
instead of the zero ^ ensemble. 

In Fig. |4] we use a theory with Np = 2 flavors of fermions, by taking the square root of the staggered fermion 
determinant. We perform reweighting from the ^ = 0 ensemble using Ri 1700 configurations. To maintain analyticity, 
in the reweighting procedure one must make sure that no cut of the complex square root function is crossed while the 
chemical potential is changed. In the complex Langevin simulations the rooting is implemented simply by multiplying 
the fermion drift terms with an appropriate factor [^ . We observe good agreement for small values n/T, similarly 
to case of the Np = 4 theory, indicating that the effect of rooting is the same in these different approaches. 


B. Reweighting from the phasequenched ensemble 

We have investigated the efficiency of reweighting from the ’phasequenched’ ensemble. On Fig. [5] we show the 
comparison of the plaquette averages as well as the Polyakov loop averages. We have used about 4000-5000 independent 
configurations at W = 4 for each /i value. One notes that the agreement is much better when compared to the 
reweighting from the /i = 0 ensemble, also for higher ^/T values (compare with Fig. [T|). Note that this comparison 
is in the deconfined phase, therefore no phase transition corresponding to the pion condensation is expected in the 
phasequenched ensemble, makeing reweighting easier. For the /3 = 5.4 value used for these plots, the complex Langevin 
simulation breaks down in the saturation region ^/T > 5 (not shown in the plots), also signaled by a large ’skirt’ 
of the distributions (meaning a slow, typically power law decay) and the disagreement of the reweighting and CLE 
simulations, most detectable in the plaquette averages. 


C. Comparisons as a function of /3 

Next we have investigated the appearance of a discrepancy of the CLE and reweighting results at smaller j3 values 
arising from a ’skirt’ of the complexified distributions [3rll31|. In Fig.|6]we compare reweighting and CLE as a function 
of the (3 parameter at fixed /r/T = 1 on an 8^ + 4 lattice. One observes that the reweighting is nicely reproduced 
by the Complex Langevin simulations as long as /3 > 5.10 — 5.15. Below this limit the distributions develop a long 
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FIG. 5: Comparison of the plaquette averages and Polyakov loops as a function of /r at a fixed /3 = 5.4 calculated with CLE 
and reweighting from the phasequenched ensemble. 
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FIG. 6: Comparison of the plaquette averages and Polyakov loops as a function of the /3 parameter at a fixed ^,/T = 1 calculated 
with CLE and reweighting from the jj, — 0 ensemble. 


skirt and CL simulations become instable, also signaled by large unitarity norm and the conjugate gradient algorithm 
(needed for the calculation of the drift terms in the CLE) failing to converge. Similar behavior is detected on the 
fermionic observables in Fig. [T] This behavior has been observed also in HDQCD simulations [13,113, where a limit 
value /3niin = 5.6 — 5.7 was seen independent of the value of Nt > 6, and /3mm was slightly smaller for Nt = 4. This 
minimal P parameter corresponds to a maximal lattice spacing amax ~ 0.2fm in HDQCD. Apparently the limiting 
/3 value is different in full QCD, but it turns out that the corresponding lattice spacing is roughly equal for Np = 4 
with am = 0.05: Umax ~ 0.2 — 0.25 fm. This breakdown is also visible on histograms of various observables. In Fig. [8] 
we show the histograms of the spatial plaquettes at various P values. One notices that the ’skirt’ of the distribution 
is indeed large at P = 5.1, where the CLE breaks down. Although a small skirt is also present at P = 5.2, it is not 
visited frequently enough to change the averages noticeably. 

A similar behavior is observed on the finer 12^ x 6 lattice, as depicted in Fig. O We used 200-300 configurations 
for the reweighting procedure on A) = 6 lattices at every P value. We observe a limiting Pmin ~ 5.15 corresponding 
to amax ~ 0.15 fm which at W = 6 allows simulations right down to the transition temperature, but not below. 

Finally we investigated Nt = 8 lattices. In Fig. [10] we show the behavior of the gauge observables, in Fig. [TTl the 
fermionic density. We used 200-300 independent configurations at each /3 value to perform the reweighting. At small 
betas the complex Langevin simulations become instable also on these lattices, which can be observed in Figs. fTOl 
and [TTl by the absence of results. One observes that the CLE breaks down above the lattice spacing a ~ 0.15fm. 
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FIG. 7: Comparison of the chiral condensate and the fermionic density as a function of the /3 parameter at a fixed /r/T = 1 
calculated with CLE and reweighting from the fj, = 0 ensemble. 
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FIG. 8: The histograms of the spatial plaquette variable in the CLE simulation, measured at various /3 values corresponding 
to Figs m and [7] 
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calculated with CLE and reweighting from the n = 0 ensemble on 12® * 6 lattices. 
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FIG. 10: Comparison of the plaquette averages and Polyakov loops as a function of the /3 parameter at a fixed /r/T = 0.96, 
using m = 0.01 and m = 0.05 as indicated, calculated with CLE and reweighting from the fi = 0 ensemble on 16® * 8 lattices. 
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m = 0.05 as indicated, calculated with CLE and reweighting from the /r = 0 ensemble on 16® * 8 lattices. 
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V. CONCLUSIONS 

In this paper we have compared complex Langevin simulations of finite density QCD with reweighting from the 
positive ensembles of the phasequenched theory and fi = 0. 

Both methods have a limited region of parameter space where they are applicable. The complex Langevin method 
fails for too small /3 parameters, as noted earlier, but this still allows the exploration of the whole phase diagram in 
HDQCD [ 4 ^. Reweighting from zero ^ breaks down because of the overlap and sign problems around ^ fa 1 — 1.5. 
In contrast, the reweighting from the phasequenched ensemble in the deconfined phase performs better also for large 
/i, suggesting that the sign problem is not that severe. 

We observe good agreement of these two methods in the region where they are both applicable. The failure of 
both methods can be assessed independently of the comparison: the complex Langevin simulations develop ’skirted’ 
distributions as the gauge cooling loses its effectiveness, and the errors of the reweighting start to grow large signaling 
sign and overlap problems. 

An important question for the applicability of the complex Langevin method to explore the phase diagram of QCD 
is the behavior of Pmim the lattice parameter below which gauge cooling is not effective. In this study we have 
determined that using W = 4, W = 6 and Nt = 8 lattices (with pion mass m-r^lTc ~ 2.2 — 2.4) this breakdown 
prevents the exploration of the deconfinement transition and the location of a possible critical point. 
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